function [TEC] = TECGauss(g1, g2, perigee, data, tol)
%
% Numerical Integration with Gauss Algorithm
%
%DESCRIPTION:
%This function implements the numerical integration of the TEC with Slat
%path using Gauss Algorithm. [2.5.8.2.8]
%
%PROTOTYPE:
% [TEC] = TECGauss(g1, g2, perigee, data, tol)
%
%--------------------------------------------------------------------------
% INPUTS:
%   g1         [1x1]       1st Integration Point     [km]
%   g2         [1x1]       2nd Integration Point     [km]
%   perigee    [---]       Ray-Perigee Data          [strc] (see NOTES)
%   data       [---]       Problem Data              [strc] (see NOTES)       
%   tol        [1x1]       Tolerance for Int. Accur. [-]    (optional)
%--------------------------------------------------------------------------
% OUTPUTS:
%   TEC        [1x1]       Total Electron Content    [TECU]
%--------------------------------------------------------------------------
%
%NOTES:
% - The input "g" is the "Distance from the Ray-Perigee of the Integration
%   Point"
% - For the tollerance (accuracy) "tol" the adviced value is:
%       between s1 and sa -> tol = 0.001
%       between sa and sb and sb and s2 -> tol = 0.01
% - The input "perigee" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       perigee.rp      = Radius           [km]
%       perigee.latp    = Latitude         [deg]
%       perigee.lonp    = Longitude        [deg]
%       perigee.sinlatp = Sine of Latitude [-]
%       perigee.coslatp = Cosine of Latit. [-]
%       perigee.sinsigp = Sine of Zenith   [-]
%       perigee.cossigp = Cosine of Zenith [-]
% - The input "data" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       data.a0      = 1st Az Coeff   [-]
%       data.a1      = 2nd Az Coeff   [-]
%       data.a2      = 3rd Az Coeff   [-]
%       data.mth     = Month          [month]
%       data.UT      = Universal Time [hours]
%       data.stModip = MODIP Table    [-] (modipneqg_wrapped.asc)
%       data.F2      = F2 Table       [-] (ccir21.asc)
%       data.Fm3     = Fm3 Table      [-] (ccir21.asc)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Accuracy Input Check
if nargin == 3
    tol = 1e-2;
end

%Nr. of Discretization Points (Initial Guess)
n = 8;

%Integration
j = 0;          %First Run Index
res = 1;        %Residual for 1st run
while (res > tol)
    if j == 1
        %Doulbling Number of Points (from 2nd run on)
        n = 2*n;
        %Definition of GN1 (past-iterate)
        GN1 = GN2;
    end
    %Integration Intervals
    Dn = (g2 - g1) / n;
    g  = 0.5773502691896 * Dn;
    y  = g1 + (Dn - g)/2;
    
    %Summatory Cycles
    aux = 0;
    for i = 1 : (n-1)
        %Determination of Altitude Case (1st Node)
        s1 = y + i*Dn;
        [N1] = ElecDens(s1, perigee, data);
        
        %Determination of Altitude Case (2nd Node)
        s2 = y + i*Dn + g;
        [N2] = ElecDens(s2, perigee, data);
        
        %Computation of the i-th Summatory Iteration
        aux = aux + (N1 + N2);
    end
    GN2 = Dn/2 * aux;
    if j == 0    %1t run
        res = 1;
    else         %From 2nd run on
        res = abs(GN1 - GN2)/abs(GN1);
    end
    %Define that 1st run is over
    j = 1;
end

%Total Electron Content
TEC = (GN2 + (GN2 - GN1)/15) * 1e-13;

end